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Abstract 

In  this  thesis,  the  problem  of  computing  the  cumulative  distribution  function 
(cdf)  of  the  random  time  required  for  a  system  to  first  reach  a  specified  reward 
threshold  when  the  rate  at  which  the  reward  accrues  is  controlled  by  a  continuous¬ 
time  stochastic  process  is  considered.  This  random  time  is  a  type  of  first  passage 
time  for  the  cumulative  reward  process.  The  major  contribution  of  this  work  is  a 
simplified,  analytical  expression  for  the  Laplace-Sticltjes  transform  of  the  cdf  in  one 
dimension  rather  than  two.  The  result  is  obtained  using  two  techniques:  i)  by  con¬ 
verting  an  existing  partial  differential  equation  to  an  ordinary  differential  equation 
with  a  known  solution,  and  ii)  by  inverting  an  existing  two-dimensional  result  with 
respect  to  one  of  the  dimensions.  The  results  are  applied  to  a  variety  of  real-world 
operational  problems  using  one-dimensional  numerical  Laplace  inversion  techniques 
and  compared  to  solutions  obtained  from  numerical  inversion  of  a  two-dimensional 
transform,  as  well  as  those  from  Monte-Carlo  simulation.  Inverting  one-dimensional 
transforms  is  computationally  more  expedient  than  inverting  two-dimensional  trans¬ 
forms,  particularly  as  the  number  of  states  in  the  governing  Markov  process  increases. 
The  numerical  results  demonstrate  the  accuracy  with  which  the  one-dimensional  re¬ 
sult  approximates  the  first  passage  time  probabilities  in  a  comparatively  negligible 
amount  of  the  time. 


IX 


TRANSIENT  ANALYSIS  AND  APPLICATIONS 
OF  MARKOV  REWARD  PROCESSES 


1.  Introduction 

1 . 1  Background, 

In  this  thesis,  the  transient  analysis  of  Markov  reward  processes  (MRPs)  will  be 
considered.  Markov  reward  processes  are  used  to  model  a  wide  variety  of  real-world 
systems.  Some  of  these  include  multiprocessor  computer  systems,  transportation 
systems,  and  manufacturing  systems.  By  modelling  such  systems  as  MRPs,  ana¬ 
lysts  are  able  to  study  performance  measures  which  would  be  difficult  to  evaluate 
otherwise.  These  performance  measures  might  be  completion  times  of  a  job,  total 
cost  associated  with  a  product,  or  the  time  required  to  transport  a  product  from  its 
origin  to  its  ultimate  destination.  Given  the  competitive  environment  that  exists  in 
today’s  world,  the  ability  to  study  and  understand  these  performance  measures  can 
be  critical  to  an  organization’s  success.  For  these  reasons,  Markov  reward  processes 
have  been  studied  extensively  in  recent  years. 

There  are  several  different  strategies  that  can  be  used  to  analyze  a  Markov 
reward  process.  First,  one  could  conduct  a  steady-state  analysis  of  the  process 
in  which  the  behavior  of  the  system  is  considered  as  time  tends  to  infinity.  In 
other  situations,  it  may  be  appropriate  to  conduct  a  transient  analysis  of  the  system 
wherein  the  cumulative  value  of  the  reward  process  is  sought  for  some  finite  time 
value. 
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1.2  Problem  Definition  and  Methodology 


The  primary  objective  of  this  thesis  is  to  investigate  the  probability  distribution 
of  the  accumulated  reward  in  a  Markov  reward  process.  The  accumulated  reward  is 
directly  influenced  by  a  random  process  that  can  be  modelled  as  a  continuous-time 
Markov  chain  (CTMC),  { Z{t )  :  t  >  0}.  The  state  of  Z(t)  controls  the  rate  at  which 
the  reward  linearly  accumulates.  It  will  be  assumed  that  the  rates  are  all  positive, 
and  therefore  the  accumulated  reward  up  to  time  t  is  a  nondecreasing  stochastic 
process.  The  properties  of  this  process  will  be  discussed  in  Chapter  3.  Finding  the 
distribution  of  the  reward  also  provides  the  distribution  of  T(x),  the  random  time 
to  accumulate  reward  level  x. 

Clearly,  a  study  of  T(x)  will  require  transient  analysis  of  the  system.  It  makes 
little  sense  to  study  the  system  as  t  — >  oo  when  the  interest  is  specifically  in  the 
system  from  time  0  to  T(x).  To  obtain  numerical  results  for  such  a  problem,  Laplace 
and  Laplace-Stieltjes  transforms  will  be  utilized.  In  some  cases,  once  the  solutions  are 
found  in  the  transform  space,  the  Laplace  transform  inversion  tables  can  be  utilized 
to  provide  an  exact  solution.  However,  it  is  more  likely  that  numerical  inversion  of 
the  transform  will  be  needed  to  provide  an  approximate  answer  to  the  problem. 

The  existing  literature  reveals  that  these  types  of  problems  are  generally  solved 
in  two-dimensional  transform  space.  However,  two-dimensional  inversion  is  not  al¬ 
ways  an  easy  (or  efficient)  task.  The  contribution  of  this  thesis  is  to  reduce  the 
dimensionality  of  these  problems,  so  that  solutions  can  be  obtained  in  one  dimen¬ 
sion  as  opposed  to  two.  Two  methods  will  be  developed  to  find  these  one-dimensional 
solutions,  and  numerical  examples  will  be  provided  demonstrating  the  utility  of  the 
solution  on  a  variety  of  operational  problems. 


1-2 


1 . 3  Thesis  Outline 


The  remainder  of  this  thesis  is  organized  in  the  following  manner.  Chapter  2 
reviews  the  previous  literature  on  the  topic.  Not  only  does  this  review  address  many 
of  the  known  fundamental  results  for  reward  processes,  but  it  also  gives  examples  of 
their  usefulness  through  some  real-world  applications.  Chapter  3  provides  a  formal 
mathematical  description  of  the  problem,  derives  both  one-  and  two-dimensional 
results,  and  discusses  numerical  inversion  of  Laplace  transforms.  In  Chapter  4,  nu¬ 
merical  examples  will  be  demonstrated,  providing  a  means  to  compare  the  one-  and 
two-dimensional  methods.  Finally,  Chapter  5  provides  conclusions,  recommenda¬ 
tions,  and  future  work. 
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2.  Review  of  the  Literature 


In  this  chapter,  a  review  of  some  relevant  literature  on  Markov  reward  models 
is  presented.  The  first  section  introduces  Markov  reward  models,  and  discuss  the 
advantages  of  Markov  modelling  over  other  alternatives.  Section  2.2  briefly  examines 
steady-state  analysis  of  Markov  reward  processes  in  the  literature.  Finally,  in  Section 
2.3,  the  transient  analysis  of  Markov  reward  processes  is  considered.  In  addition,  this 
section  provides  examples  from  the  literature  of  specific  systems  modelled  as  Markov 
reward  processes,  and  discusses  performance  measures  and  the  methods  with  which 
they  are  obtained. 

2.1  Markov  Reward  Processes 

Markov  reward  models,  also  referred  to  as  Markov  reward  processes  (MRP), 
allow  analysts  to  accurately  model  systems  that  evolve  stochastically  over  time. 
A  Markov  reward  model  consists  of  two  elements:  a  structure-state  process  and 
a  reward  structure.  The  structure-state  process  is  assumed  to  be  an  irreducible, 
positive  recurrent  continuous-time  Markov  chain  (CTMC)  with  a  specified  state 
space,  S.  The  reward  structure  consists  of  the  rates  at  which  reward  (or  cost)  is 
accrued  when  the  structure-state  process  occupies  each  state  in  S.  Thus,  the  reward 
rate  depends  explicitly  on  the  state  of  the  system. 

One  example  of  such  a  system  is  presented  by  Kulkarni  [13].  Consider  a  ma¬ 
chine  that  can  be  in  one  of  two  states:  it  is  either  functioning  or  it  is  not  functioning. 
If  it  is  functioning,  it  fails  after  a  random  amount  of  time  which  is  exponentially  dis¬ 
tributed  with  rate  fi.  If  it  is  not  functioning,  it  is  repaired  after  a  random  amount 
of  time  which  is  exponentially  distributed  with  rate  A.  Next,  define  a  random  vari¬ 
able  X(t),  such  that  if  X{t)  =  0,  the  machine  is  not  functioning  at  time  t  and  if 
X(t)  =  1,  the  machine  is  functioning  at  time  t.  Therefore,  the  state  of  the  machine 
is  described  by  the  continuous-time  Markov  chain  (X(f)  :  t  >  0},  which  is  called 
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the  structure-state  process.  Now  suppose  that  the  machine  produces  profit  at  rate 
r  dollars  per  unit  time  when  it  is  functioning.  But  when  the  machine  is  down,  the 
repairs  cost  c  dollars  per  unit  time.  This  constitutes  the  reward  structure  in  this 
model.  When  X(t)  =  0,  the  rate  is  — c,  and  when  X(t)  =  1,  the  rate  is  r. 

Kulkarni  [13]  gives  another  example  which  involves  a  machine  shop  that  has 
two  machines  that  are  independent,  identical,  and  have  the  same  failure  and  repair 
rates  as  the  machine  in  the  previous  example.  This  time,  let  X(t)  be  the  number 
of  working  machines  at  time  t.  It  is  easily  shown  that  the  structure-state  process, 
(X(t)  :  t  >  0},  is  a  continuous-time  Markov  chain.  As  one  would  expect,  the  rate  at 
which  the  machine  shop  produces  work  at  a  specified  time,  t,  is  completely  dependent 
on  the  number  of  functioning  machines  at  time  t.  Therefore,  there  exist  different 
reward  rates  when  there  are  zero,  one,  or  two  machines  operating.  This  is  another 
example  of  a  Markov  reward  process. 

Once  a  system  has  been  modelled  analytically  as  a  Markov  reward  process, 
the  model  may  be  analyzed  to  determine  various  performance  measures,  such  as 
the  expected  time  to  complete  a  job,  or  the  expected  time-averaged  cost.  Although 
analytical  modelling  is  not  the  only  technique  used  to  make  predictions  about  system 
performance,  it  is  often  the  preferred  technique.  Besides  analytical  modelling,  there 
are  two  other  general  methods  for  making  predictions  about  the  performance  of  a 
system  [18].  The  first  technique  involves  actual  observation  and  measurement  of 
performance  metrics.  Clearly,  repeated  measurements  over  an  extended  period  of 
time  on  the  actual  system  being  studied  would  help  analysts  assess  the  behavior  of 
a  given  performance  measure.  However,  there  are  two  major  disadvantages  to  this 
technique.  Unfortunately,  it  is  often  quite  expensive  to  collect  the  data.  Moreover, 
actual  measurements  cannot  be  collected  if  the  system  does  not  exist.  As  an  example, 
if  the  study  is  performed  to  select  the  most  efficient  design  of  an  assembly  line,  this 
technique  cannot  be  utilized,  since  the  assembly  line  does  not  exist. 
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Another  alternative  for  making  predictions  is  computer  simulation.  There  are 
many  computer  software  packages  available  to  simulate  a  wide  variety  of  systems 
ranging  from  the  very  simple  to  the  very  complex.  But  as  Reibman,  et  al.  [18]  point 
out,  simulation  can  also  be  expensive  and  may  require  considerable  time  to  produce 
statistically  significant  results. 

Because  of  the  noted  disadvantages  of  the  other  two  techniques,  analytical 
modelling  is  generally  used  whenever  possible.  There  are  several  types  of  analyti¬ 
cal  models  found  in  the  literature,  including  combinatorial  models,  Markov  chains, 
Markov  renewal  processes,  and  Markov  reward  processes.  Markov  models  are  able 
to  capture  complex  attributes  of  a  system’s  behavior  [18],  and  will  therefore  be  the 
focus  of  this  literature  review. 

One  of  the  major  advantages  of  Markov  reward  processes  over  other  types  of 
analytical  models  is  that  they  allow  the  analyst  to  combine  traditional  performance 
with  reliability.  In  the  past,  models  have  predicted  system  performance  under  the 
assumption  that  systems  operate  failure-free  [16].  Obviously,  real-world  systems 
cannot  operate  indefinitely  without  experiencing  failures.  These  failures  (and  the 
associated  repairs)  are  often  modelled  using  reliability  analysis.  However,  in  recent 
years,  analysts  have  become  increasingly  aware  of  the  fact  that  the  separation  of 
performance  and  reliability  models  is  no  longer  adequate  [16].  This  has  led  to  an 
increased  interest  in  applying  Markov  reward  modelling  to  systems  in  which  the 
structure-states  represent  various  levels  of  degradation.  This  application  is  particu¬ 
larly  prevalent  in  the  area  of  computer  systems  [19]. 

It  is  evident  from  the  literature  that  Markov  reward  processes  can  be  an  ex¬ 
tremely  useful  tool.  However,  in  building  a  Markov  model,  the  analyst  should  con¬ 
sider  five  major  aspects  to  maximize  its  usefulness  [18]: 
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•  Make  the  model  easy  to  use.  If  the  model  is  too  complicated,  the  advantages 
of  Markov  modelling  may  be  outweighed  by  the  disadvantage  of  its  difficulty 
to  use. 

•  Minimize  computation  costs.  A  model  is  much  more  useful  if  does  not  require 
tremendous  amounts  of  computation  time,  or  computer  memory. 

•  Error  analysis.  An  analyst  should  understand  the  sources  of  error,  and  seek 
to  minimize  that  error.  Some  examples  of  error  sources  are  modelling  assump¬ 
tions,  estimated  parameter  values,  and  numerical  approximations. 

•  Process  improvement.  In  developing  a  model,  an  analyst  learns  a  great  deal 
about  the  process,  and  can  often  identify  areas  in  which  the  process  could  be 
improved. 

•  Identify  and  compute  desired  measures.  This  is  generally  the  purpose  of  de¬ 
veloping  the  model  in  the  first  place,  and  is  therefore  the  primary  goal. 

Once  the  model  has  been  built,  two  types  of  analyses  can  be  performed:  steady- 
state  and  transient.  In  a  steady-state  analysis,  one  considers  the  behavior  of  the 
system  as  time  tends  toward  infinity.  Predictions  of  system  performance  are  subse¬ 
quently  based  on  this  limiting  behavior.  On  the  other  hand,  a  transient  analysis  is 
concerned  with  the  system  behavior  in  the  short  run.  Therefore,  predictions  must 
be  based  on  the  behavior  of  the  system  after  a  finite  amount  of  time.  The  next  sec¬ 
tion  will  discuss  some  of  the  known  steady-state  results  for  Markov  reward  processes 
found  in  the  literature. 

2.2  Steady-State  Analysis  of  Markov  Reward  Processes 

An  extremely  important  measure  in  steady-state  analysis  is  the  steady-state 
probability  vector,  p,  which  describes  the  long-run  chance  of  finding  the  process  in  a 
given  state.  To  find  p  for  a  CTMC  (X(f)  :  t  >  0},  it  is  necessary  to  introduce  some 
notation.  First  define  qtJ  as  the  rate  at  which  the  CTMC  transitions  from  state  i  to 
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state  j  where  i  ^  j  because  a  “transition”  is  defined  as  moving  into  a  different  state. 
Further,  define 

Qii  ^  ^  Qij  i 
j+i 


and  define  the  square  matrix 


Q  —  [Qij]- 


The  matrix  Q  is  called  the  infinitesimal  generator  matrix  of  (A"(f)  :  t  >  0}.  Then 
let 

Pj  =  ,lim  pix(t)  =j} 

t—>  oo 

denote  the  long-run  probability  that  the  CTMC  is  in  state  j.  It  is  well-known  [13], 
[19],  that  for  an  irreducible,  positive  recurrent  CTMC,  the  steady-state  distribution 
is  given  by  the  the  unique  solution  to 


pQ  =  0  -  1,  (2.1) 

jes 

where  S  is  the  state  space. 

Analysts  can  use  the  steady-state  distribution  of  Equation  (2.1)  to  make  pre¬ 
dictions  about  a  system’s  performance  in  the  asymptotic  region.  For  example,  let 
rx(t )  denote  the  reward  rate  at  time  t,  and  let 

W(t)  =  1  fx(T)dT 

be  the  time-averaged  accumulated  reward  in  the  system  for  t  >  0.  Smith,  et  al.  [19] 
show  that  in  the  limit  as  t  tends  to  infinity,  the  expected  values  of  rX(t)  and  W (t) 
are  equivalent,  and  are  expressed  as 


lim  E[rX(t)]  =  bm  E[W(t)]  =  Vr®. 

t — >oo  t—>oo  z ' 

i 
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Although  much  information  can  be  gleaned  from  steady-state  analysis,  it  is 
often  more  useful  to  conduct  a  transient  analysis.  Obviously,  in  real  applications, 
many  of  the  measures  that  analysts  seek  involve  some  measure  up  to  a  time  t.  For 
example,  a  first  passage  time  is  the  random  time  required  for  a  system  to  first  accu¬ 
mulate  a  specified  reward  level.  Clearly,  it  makes  no  sense  to  use  limiting  behavior 
in  such  a  situation. 

2.3  Transient  Analysis  of  Markov  Reward  Processes 

A  review  of  the  literature  on  transient  analysis  of  Markov  reward  processes 
shows  that  MRPs  have  wide  applicability  in  the  real  world.  In  this  section,  several 
examples  from  the  literature  will  be  provided  to  show  not  only  the  variety  of  types 
of  systems  that  can  be  modelled  with  a  Markov  reward  process,  but  also  to  show  the 
variety  of  analysis  approaches  that  exist. 

The  first  example  comes  from  [19],  in  which  Smith,  et  al,  present  a  case  study 
in  which  they  model  a  multiprocessor  system  with  16  processors,  16  memories,  and  a 
crossbar  switch.  The  “structure-state”  of  the  system  at  any  given  time  is  represented 
by  the  triple  ( i,j,k )  indicating  the  number  of  operational  processors,  memories, 
and  switches,  respectively.  For  the  system  to  be  functioning,  the  switch  must  be 
operational,  and  a  certain  number  of  processors  and  memories  must  be  functional. 
In  this  paper,  they  specify  that  there  must  be  four  processors  and  four  memories  in 
operation.  Therefore,  this  system  has  169  functioning  states  and  1  failed  state  for  a 
total  of  170  states.  The  authors  also  point  out  that  by  altering  the  types  of  failure 
transitions,  and  the  types  of  repairs,  the  number  of  states  can  grow  to  as  much  as 
365. 

Having  established  the  structure-state  process,  the  authors  consider  the  reward 
structure  for  this  example.  They  present  three  possibilities.  The  first  is  a  simple 
availability-based  structure  in  which  all  operational  states  are  assigned  a  reward  rate 


2-6 


of  1,  and  the  failure  state  is  assigned  a  reward  rate  of  0.  A  second,  more  accurate, 
structure  is  a  capacity-based  structure  in  which  each  operational  structure-state 
(i,j,  1)  is  assigned  reward  rate  min{i,j},  and  the  failure  state  is  assigned  rate  0. 
The  third  possible  structure,  developed  by  Bhandarkar  [4],  dictates  that  the  reward 
rate  in  an  operational  state  (i,j,  1)  is  rtJj  =  m(l  —  (1  —  1  /m)L),  where  L  =  min  {7,  j} 
and  m  =  max{i,j}.  Again,  the  failure  state  is  assigned  rate  0. 

Next,  the  authors  examine  performability  results,  based  on  the  170  state 
structure-state  process,  and  the  contention-based  reward  structure.  With  the  neces¬ 
sary  elements  of  the  Markov  reward  model  defined,  they  identify  and  compute  four 
different  measures.  The  first  is  the  computation  of  availability,  E[X(t)\,  where  X{t) 
is  defined  as  the  reward  rate  at  time  t,  which  answers  the  question  “What  is  the 
expected  performance  of  the  system  at  time  £?”  Next,  they  identify  the  expected 
time-averaged  accumulated  reward  over  the  interval  (0,t),  which  answers  the  ques¬ 
tion  “What  is  the  time-averaged  performance  of  the  system  over  the  interval  (0,t)?” 
The  third  measure  they  identify  is  the  likelihood  of  completing  a  given  amount  of 
work  in  a  specified  time  interval,  yc(x,t),  which  answers  the  question  “What  is 
the  probability  that  x  units  of  work  are  completed  by  time  £?”  The  final  measure 
is  Wc(x,t),  which  answers  the  question  “What  is  the  probability  that  the  reward 
accumulated  in  the  interval  (0,  t )  is  at  least  xtE 

Because  yc(x,t )  is  a  first  passage  time,  it  is  of  particular  interest  to  this  thesis. 
In  [19],  the  authors  derive  an  expression  for  yc(x,t )  in  two-dimensional  transform 
space  which  requires  numerical  inversion  to  find  values  for  the  cdf  of  the  first  pas¬ 
sage  time  for  specified  values  of  t.  This  transform  solution  is  similar  to  the  two- 
dimensional  transform  solution  reviewed  in  Chapter  3  of  this  thesis. 

Another  example  of  a  real-world  computer  system  modelled  as  a  Markov  reward 
model  is  given  by  Kulkarni,  et  al.  [12],  The  authors  present  a  general  model  of  the 
completion  time  of  a  single  job  on  a  computer  system  whose  state  changes  according 
to  a  Markov  process.  Unlike  the  previous  model,  when  the  state  of  the  system 
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changes,  the  service  is  preempted.  The  service  is  then  resumed  or  restarted  in  the 
new  state  at  a  rate  that  could  possibly  be  different  than  the  original  rate.  Thus,  this 
model  incorporates  three  types  of  service  disciplines:  preemptive-resume  (PRS),  in 
which  the  job  is  restarted  from  the  point  where  it  was  preempted;  preemptive-repeat- 
identical  (PRI),  in  which  the  job  is  restarted  from  the  beginning;  and  preemptive- 
repeat-different  (PRD),  in  which  the  job  is  restarted  with  a  new  work  requirement. 

In  their  model,  they  define  B  >  0  as  the  amount  of  work  required  to  complete 
a  job,  and  { Z{t )  :  t  >  0}  as  the  continuous-time  stochastic  process  defined  on  state 
space  S  =  {1,  2, ...}.  The  structure-state  of  the  system  is  described  by  this  stochastic 
process.  Each  state  is  also  identified  as  either  PRS,  PRI,  or  PRD.  The  reward  rate 
of  the  system  is  such  that  when  the  system  is  in  state  i,  the  work  rate  is  ry  >  0. 
Finally,  T  is  defined  as  the  time  needed  to  complete  a  job  with  requirement  B.  It  is 
the  cumulative  distribution  function  of  the  random  variable  T  that  the  authors  seek. 

To  End  this  cdf,  Kulkarni,  et  al.  [12]  take  a  unique  approach  which  they  call 
“progressive  aggregation.”  This  method  is  described  as  a  12-step  process  involving 
a  series  of  transforms  and  inversions  to  ultimately  arrive  at  an  expression  for  the 
Laplace-Stieltjes  transform  of  F,  which  is  the  cdf  of  T.  This  expression  is  a  sum  of 
i  weighted  Laplace-Stieltjes  transforms,  where  i  is  the  number  of  functional  states 
in  the  system.  The  authors  demonstrate  their  technique  on  a  system  with  four 
structure-states  {0, 1,  2,  3},  where  state  0  is  a  failure  state,  state  1  is  PRS,  state  2  is 
PRI,  and  state  3  is  PRD.  The  structure-state  process  is  a  CTMC  with  A ij  denoting 
the  rate  of  transition  from  state  i  to  state  j,  i  ^  j.  The  reward  rates  for  states  1,  2, 
and  3  are  ry,  r2,  and  r3,  respectively  (r0  =  0).  Following  their  12-step  process,  they 
finally  arrive  at  the  solution 

3 

F«  =  £P{Z(0)  =.}£«, 

i= 1 

where  F,:(s)  is  the  Laplace-Stieltjes  transform  of  F)(f)  =  P{T  <  f|Z(0)  =  i}. 
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Although  the  authors  arrive  at  a  solution  in  one- dimensional  transform  space, 
computation  could  be  difficult  and  extremely  time-consuming,  particularly  as  the 
number  of  states  increases.  This  problem  differs  from  the  one  addressed  in  this 
thesis;  however,  it  demonstrates  that  it  is  possible  to  accurately  model  very  complex 
systems  using  a  Markov  reward  model. 

In  the  literature,  there  are  other  examples  of  computer  systems  being  modelled 
by  Markov  reward  models,  with  many  researchers  finding  a  solution  to  the  cdf  of 
the  first  passage  time  as  a  two-dimensional  Laplace-Stieltjes  transform  (as  seen  in 
Chapter  3).  Occasionally,  one  finds  a  unique  approach  such  as  the  previous  example. 
Another  exception  is  the  model  developed  by  Nicola,  et.al.  [16]  in  which  the  cdf 
is  expressed  as  a  functional  equation  involving  a  one-dimensional  Laplace-Stieltjes 
transform.  Solving  this  equation  can  be  difficult,  but  it  provides  an  alternative  to 
the  two-dimensional  approach. 

Clearly,  computer  systems  are  excellent  examples  of  these  types  of  models,  but 
Markov  reward  models  are  certainly  not  limited  to  computer  systems.  Kharoufeh 
[8]  provides  a  different  type  of  example.  In  his  dissertation,  he  presents  a  variant 
of  a  Markov  reward  model  to  analyze  a  problem  from  vehicular  traffic  flow  theory. 
Consider  a  vehicle  travelling  along  a  roadway  segment.  The  various  states  of  the 
structure-state  process  are  represented  by  the  environmental  conditions  to  which 
the  vehicle  is  subjected  during  its  sojourn.  For  example,  one  state  might  coincide 
with  rainy  conditions,  and  another  might  indicate  clear  conditions.  The  state  of  the 
environment  at  any  given  time  is  governed  by  a  continuous-time  Markov  chain.  The 
reward  rates  in  this  system  are  the  velocities  at  which  the  vehicle  may  travel,  and 
are  determined  by  the  state  of  the  system  such  that  when  the  environmental  process 
is  in  state  i,  the  vehicle  travels  at  velocity  Vt.  Finally,  the  accumulated  reward  is  the 
total  displacement  of  the  vehicle  up  to  a  fixed  point  in  time.  Kharoufeh  [8]  seeks  to 
find  several  measures  in  this  system.  The  primary  measure  of  interest  to  this  thesis 
is  the  unconditional  cumulative  distribution  function  of  the  first  passage  time  T(x), 
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where  T(x)  is  the  random  time  required  to  first  travel  a  distance  of  x.  Kharoufch’s 
solution,  which  is  examined  in  detail  in  Chapter  3,  is  expressed  as  a  two-dimensional 
transform  of  the  cdf,  and  must  be  inverted  numerically  to  obtain  results.  In  addition 
to  the  cdf  of  the  first  passage  time,  the  author  derived  the  transient  and  asymptotic 
moments  of  the  first  passage  time. 

Markov  reward  models  can  be  advantageous  to  other  alternatives  for  an  as¬ 
sortment  of  reasons.  They  are  generally  less  costly  in  terms  of  time  and  money 
than  the  alternative  methods.  Also,  the  reward  rates  in  a  Markov  reward  pro¬ 
cess  provide  modelers  with  a  technique  to  easily  incorporate  reliability  into  their 
models,  making  Markov  reward  processes  increasingly  more  popular.  This  chapter 
presented  some  specific  examples  found  in  the  literature,  as  well  as  the  pertinent 
performance  measures.  In  these  examples,  it  was  seen  that  multiple  techniques  have 
been  employed  to  find  the  cdf  of  various  first  passage  times  for  the  reward  pro¬ 
cess.  Some  involved  solutions  that  required  numerical  inversion  of  two-dimensional 
Laplace  transforms,  and  others  required  extensive  calculations  to  arrive  at  a  solution 
requiring  one-dimensional  Laplace  inversion.  However,  the  existing  literature  does 
not  contain  a  simplified,  analytical  expression  for  the  Laplace-Stieltjes  transform  of 
the  cdf  in  one  dimension.  The  goal  of  this  thesis  is  to  provide  a  different  approach 
to  develop  an  expression  for  the  first  passage  time  distribution  which  can  be  found 
easily  and  with  very  little  computation  time. 
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3.  Formal  Mathematical  Model 

The  primary  objective  of  this  chapter  is  to  present  a  detailed  description  of 
Markov  reward  processes  and  to  provide  an  explicit  result  for  the  first  passage  time 
distribution.  Additionally,  numerical  methods  for  computing  approximations  to  the 
distributions  will  be  explored. 

3.1  Markov  Reward  Processes 

Let  {Z(t)  :  t  >  0}  be  a  continuous-time  Markov  chain  (CTMC)  having  finite 
state  space  S  =  {1,2 for  some  value  K  G  N,  the  set  of  positive  integers. 
Furthermore,  let  Q  =  [qy,-],  i,j  G  S,  denote  the  infinitesimal  generator  matrix  for 
{Z{t)  :  t  >  0}  such  that  qy,-  denotes  the  rate  at  which  the  process  transitions  from 
state  i  G  S  to  j  G  S,  i  ^  j  and  qlt  =  -  Qij- 

The  evolution  of  a  Markov  reward  process  can  be  described  as  follows.  At  time 
0,  the  stochastic  process,  {Z(t)  :  t  >  0},  has  initial  distribution  z0  and  the  total 
accumulated  reward  (or  cost)  is  zero.  Whenever  Z(t)  assumes  a  value  of  i  G  S',  the 
system  accumulates  reward  at  rate  rt  >  0.  The  total  accumulated  reward  up  to  time 
t  is  represented  by  the  random  variable  R(t).  Therefore,  the  process  {R(t)  :  t  >  0} 
is  a  continuous-time  stochastic  process  on  the  continuous  space  [0,  oo).  The  reward 
could  be  a  distance  travelled,  monetary  income,  or  even  a  cost  such  as  accumulated 
damage  to  a  component.  For  the  purposes  of  this  thesis  it  will  be  assumed  that  the 
reward  is  non-decreasing.  Because  the  state  space  is  finite,  there  is  also  a  finite  set 
of  reward  rates,  1Z  =  {ri,  r2, ...,  tk}-  The  accumulated  reward  in  the  process  can 
then  be  expressed  as 

R(t)  =  R{ 0)  +  [  rz(u)du , 

Jo 

where  R(0)  is  the  reward  level  at  time  0,  which  will  be  assumed  to  be  0  throughout. 
As  an  example  of  a  Markov  reward  process,  Figure  3.1  demonstrates  the  relation- 
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ships  between  the  various  random  variables  in  a  sample  three-state  system.  In  this 
example,  rt  =  3  —  i. 

Now  let  Pij(t )  =  P{Z(t)  =  j\Z(0)  =  i}  be  the  probability  that  the  CTMC  is 
in  state  j  at  time  t,  given  that  the  CTMC  was  in  state  i  at  time  0,  where  i,j  G  S. 
Further,  define  P(t)  =  [pij(t)\,  the  transition  probability  matrix,  which  is  described 
by  the  system  of  forward  Kolmogorov  equations 

<^L  =  P(t)Q ,  Po  =  P{  0).  (3.1) 

The  well-known  solution  to  the  initial  value  problem  in  (3.1)  is 

Pit)  =  P0exp  (Qt),  (3.2) 

where  exp  (Qt)  is  defined  by  the  infinite  series 

exp  (Qt)  =  /  +  gt+^r  +  ^r  +  ...  (3.3) 

that  is  convergent  for  every  choice  of  Q  and  t. 

Finally,  let  R(t)  be  the  total  reward  up  to  time  t,  and  define 

T(x)  =  inf {t  :  FL(t)  >  a:}, 

the  random  time  it  takes  to  first  achieve  a  total  reward  of  x  G  M+.  The  random  time, 
T(x ),  is  often  referred  to  as  the  “first  passage  time”  for  the  process  {R(t)  :  t  >  0}. 
The  main  objective  of  this  thesis  is  to  explore  numerical  methods  for  the  evaluation 
of  the  cumulative  distribution  function  of  T(x ). 
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Figure  3.1  Markov  reward  process:  three-state  CTMC. 
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3.2  Distribution  of  the  First  Passage  Time:  Two-Dimensional  Ap¬ 
proach 

In  this  section,  basic  probability  principles  will  be  used  to  obtain  the  distribu¬ 
tion  of  T(x)  in  two-dimensional  transform  space.  First,  it  is  necessary  to  define  the 
joint  probability  distribution 

Hfx,t )  =  P[R(t)  <  x,Z(t )  =  i] 


so  that  by  Bayes’  Rule, 

P[R(t)<x]  =  J2Hi(x,t). 

i£S 

However,  it  can  be  shown  that  under  the  assumption  that  R(t)  is  non-decreasing, 
the  event  (-R(t)  <  a:}  is  equivalent  to  the  event  {T(x)  >  t}.  Thus,  if  Gx{t )  is  the 
cumulative  distribution  function  (cdf)  of  T(x),  then 


Gx(t)  =  P[T{x )  <t\  =  l-Y,  *)■  (3-4) 

i£S 


Kharoufeh  and  Gautam  [9]  derived  the  distribution  of  T(x)  by  showing  that 
(. x ,  t)  satisfies  the  partial  differential  equation  (PDE) 


d  Hi(x,t) 

dt 


+ 


d  Hj(x,t) 

dx 


r<  = 


J ~2qjiHi(x,t ), 

ies 


i  e  S 


(3.5) 


where  rt  E  7 Z.  Additionally,  set  the  initial  condition  Hfx,  0)  =  P[Z(0)  =  i\.  This 
equation  can  then  be  converted  to  matrix  form  by  letting  H(x,t )  =  [Hfx,t)]i£S  be 
the  row  vector  and  V  =  diag(ri,  r2,  ...,  rK).  Equation  (3.5)  then  becomes 


d H(x,  t ) 
dt 


+ 


d H(x,  t ) 
dx 


V  =  H(x,t)Q . 


(3.6) 


for  t  >  0  and  x  >  0. 


3-4 


Next,  define  zq  as  the  initial  distribution  of  Z(t).  Also,  define 


H*(x,s2)  =  /  e  S2tH(x,t)dt, 

Jo 

the  Laplace  transform  of  H(x,  t)  with  respect  to  t,  and 


H*(s1,s2)=  /  e~SlXdH*(x,s2), 


Jo 

the  Laplace-Stieltjes  transform  of  H*(x,s2)  with  respect  to  x.  Kharoufeh  and  Gau- 
tam  [9]  proved  that  the  solution  to  Equation  (3.6)  in  the  transform  domain  is 

H*(si,  s2)  —  Zq(siV  +  s2I  —  Q)  1  (3-7) 


where  si  and  s2  are  complex  variables  with  Re(s i)  >  0  and  Re(s2 )  >  0.  Equation 
(3.7)  can  be  solved  numerically  using  a  two-dimensional  inversion  algorithm.  Moor- 
thy  [15]  and  Abate,  et  al.  [3]  provide  two  such  algorithms,  which  will  be  discussed 
later  in  this  thesis. 

When  numerically  inverting  Equation  (3.7),  one  occasionally  runs  into  diffi¬ 
culties  due  to  stability  issues  relating  to  the  inversion  algorithm.  This  issue  will  be 
discussed  further  in  Section  3.4.  Additionally,  the  numerical  inversion  can  often  be 
time  consuming  in  complex  problems  in  which  the  CTMC  has  a  large  number  of 
states,  (K  >  10).  Because  of  these  problems,  it  would  be  beneficial  to  reduce  the 
problem  to  a  single  dimension  transform.  Numerical  inversion  in  one  dimension  can 
be  considerably  quicker,  and  does  not  suffer  from  the  same  instability  problems. 


3.3  Distribution  of  the  First  Passage  Time:  One- Dimensional  Ap¬ 
proach 

The  main  result  of  this  thesis  provides  a  reduction  in  the  dimensionality  of  the 
first  passage  time  distribution  of  a  Markov  reward  process  from  two  dimensions  to 
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a  single  dimension  and  is  stated  in  Theorem  3.1.  The  result  shall  be  proved  in  two 
ways.  First  via  direct  methods  and  second,  by  inverting  the  result  of  Kharoufeh  and 
Gautam  [9]  with  respect  to  x. 

Theorem  3.1  Let  Gx(s )  denote  the  Laplace- Stieltjes  transform  of  the  cumulative 
distribution  function  ofT(x).  Then, 

Gx(s )  =  z0exp(E_1[Q  -  sl]x)l,  (3.8) 


where  1  —  [1, 1,  •  •  •  ,  1]T. 

3.3.1  Proof  of  Theorem  3.1:  Method  One 

In  this  method,  the  desired  result  is  obtained  directly  from  Equation  (3.6), 


mx,t)  +  mxjlv  =  H(xt)Q 


dt 


dx 


Taking  the  Laplace  transform  of  Equation  (3.6),  with  respect  to  t  yields 


dH*(r  si 

sH*(x,  s)  -  H(x,  0)  +  - — =  H*(x,  s)Q, 

dx 


(3.9) 


since 

£  (^HQf  =  sH*(x,  s )  -  H{x,  0). 

But  Hfx,  0)  =  P{R( 0)  <  x,Z( 0)  =  i}  =  P{Z{ 0)  =  i}  and  therefore  H(x,  0)  =  zo, 
where  Zq  is  the  initial  distribution  vector  of  Z(t).  Equation  (3.9)  then  becomes 

dH*(r  si 

sH*(x,s)-z0  +  - — f^-V  =  H*(x,s)Q.  (3.10) 

dx 
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Rearranging  terms  leaves  a  linear  ordinary  differential  equation  (ODE)  system  with 
constant  coefficients, 


dH*^s)  +H*{xJs)[sI-Q}V~l  =  z0V-\  (3.11) 

For  this  ODE,  the  necessary  integrating  factor  is 

,(n=exp(/[s/-Q]^^)=exp([s/-Q]^n). 

since  [si  —  Q\V~l  is  a  constant  matrix  with  respect  to  x.  Right  multiplying  both 
sides  of  Equation  (3.11)  by  yields 

(. H*(x ,  s)  exp([s/  —  Q^V^x))'  =  z0V -1  exp([s/  —  Q]V~lx).  (3.12) 

Integrating  both  sides  of  Equation  (3.12)  with  respect  to  x  produces 

H*(x,  s )  exp([s/  —  Q]V~lx)  =  ZqV-1  (V[sI  —  Q]_1)  exp([sJ  —  Q]V~lx)  +  C,  (3.13) 

where  C  is  the  constant  row  vector  of  integration.  Simplifying  and  applying  the 
initial  condition  H*(0,s)  =  0  shows  that  C  =  —  zq\sI  —  Q]~l .  This  results  in 

H*(x,s)  =  z0[sl  -  QY1  -  zq[sI  -  QY1  exp([Q  -  sljV^x),  (3.14) 


which  can  be  shown  to  be  equivalent  to 


H*(x,  s)  =  zq[sI  —  Q]  *  —  zoexp(E  l[Q  —  sl]x)[sl  —  Q]  1.  (3.15) 


Define  the  Laplace-Sticltjes  transform  of  the  distribution  of  the  first  passage 


time  as 


£[Gx](s)  =  Ga(s)=  /  e~stdGx(t), 
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which  is 


Gx(s)  =  1  —  H*(x,  s)sl,  (3.16) 

where  1  is  a  column  vector  of  ones.  Combining  Equations  (3.14)  and  (3.16)  produces 

Gx(s)  =  1  -  (z0  -  zo  exp(U_1[<5  -  sl]x))[sl  -  <2]-1sl.  (3.17) 

This  reduces  to 

r  or1 

Gx(s)  =  1  -  (z0  -  z0exp(U_1[Q  -  sl]x))  1.  (3.18) 

Using  the  Neumann  expansion  [17]  (I  —  A)~l  =  /  +  A  +  A2  +  A3  +  •  •  • ,  applied  to 
A  =  Q/s,  where  ||A||  =  max  ]T)  | aA  <  1,  Equation  (3.18)  becomes 

*  3 

Gx(s )  =  1  -  (1  -  zo  exp(U^1[Q  -  sl}x)  1),  (3.19) 

since  Q 1  =  0  and  zqI  =  1  for  an  irreducible,  finite-state  CTMC.  This  gives  the 
desired  result 

Gx(s )  =  ^0exp(U_1[<5  —  sl]x)  1.  (3.20) 

3.3.2  Proof  of  Theorem  3.1:  Method  Two 

In  the  second  method,  the  desired  one- dimensional  result  is  obtained  from  the 
two-dimensional  result  in  Equation  (3.7). 

From  Equation  (3.16), 

Gx(s)  =  1  —  H*(x,  s)sl. 
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Taking  the  Laplace-Sticltjes  transform  of  both  sides  with  respect  to  x,  and  combining 
the  result  with  Equation  (3.7)  yields 

G(w)  =  1  —  H*(w ,  s)sl  =  1  —  zq(wV  +  si  —  Qy'ls, 

where  G(w)  represents  the  double  Laplace-Sticltjes  transform  of  G(x).  Since  Laplace 
transforms  are  easier  to  invert,  the  above  equation  is  written  in  terms  of  the  Laplace 
transform  as 

wG*(w )  =  1  —  zq(wV  +  si  —  Q)_1ls, 
because  G  is  related  to  G*  by  the  equation  [13] 

G(w)  =  wG*(w). 


Rearranging  terms  results  in 


wG*(w )  =  1  —  — 
w 


I  - 


V~\Q-sI) 


w 


-l 


V-1!  s. 


Again  using  the  geometric  expansion  (/  —  A)  1  =  I  +  A  +  A2  +  A3  +  •  ■  ■ ,  for  ||  A||  <  1, 
and  dividing  by  w  gives 


G*(w)  = 


1  z0 


w 


n=  1 


V~\Q-sI) 


w 


V~lls, 


for  ||  V  1{Q  —  s/)||  <  1.  In  order  to  facilitate  inversion,  the  above  equation  is  rewritten 


as 


G*{w)  = - z0 

w 


wz 


oo  1 

Erw-oi-A 


n= 1 


V^ls. 
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Since  the  inverse  Laplace  transform  is  uniformly  continuous,  the  transform  of  the 
infinite  sum  is  the  infinite  sum  of  the  transforms.  That  is, 

OO  1~00  f  1 

£  -  si)r—2  =^(v-h« ^  • 

_n=  1  J  n=  1  L  J 

The  terms  of  this  series  can  now  be  inverted,  resulting  in 

OO  |_^ 

<?*(«)  =  +  V-'ls. 

n= 1  '  '  _ 

After  some  algebra,  the  above  equation  can  be  written  as 

(V~1(D  —  <?  T)r)n+1 

Gx(s)  =  l-z0  xV-\Q-sI)  +  y2[- - 7  ,  (Q  —  sl)~1ls.  (3.21) 

“  (n+1)! 

71=  1  x  7 

Since  the  exponential  of  a  matrix  is  defined  as  exp(A)  =  77o  A? /j\.  Equation  (3.21) 
can  be  written  as 

Gx(s)  —  l-z0  [-1  +  exp(E_1(Q  -  sl)x)]  ( Q  -  s/)-1ls. 

Once  again,  expanding  (Q  —  sl)^1  by  the  Neumann  expansion  yields 


ffowever,  because  01  =  0  and  z0 1  =  1,  the  result  is 

Gx(s )  =  z0exTp(V~1(Q  —  sl)x)  1.  (3.22) 


Gx(s )  —  l  +  ^o  I  +  exp(E  1(Q  —  s/)x)]  f  I  H - 1  —  H  — 

V  s  s  s 


3.4  Numerical  Inversion 

Once  a  solution  is  obtained  in  the  transform  domain,  the  problem  of  inverting 
that  solution  remains.  In  certain  cases,  an  exact  expression  for  the  inverse  transform 
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can  be  obtained.  For  example,  if  the  transform  solution  is  a  vector  of  rational  func¬ 
tions  in  all  the  complex  variables,  partial  fraction  decomposition  can  be  employed 
to  obtain  a  form  which  facilitates  analytic  inversion  via  the  inverse  Laplace  tables. 
However,  this  can  be  a  difficult  task.  Therefore,  obtaining  approximations  by  nu¬ 
merical  inversion  is  generally  preferred.  At  this  point,  there  will  be  a  brief  review  of 
numerical  inversion  of  Laplace  transforms  in  both  one  and  two  dimensions. 


3-4-1  Numerical  Inversion  in  One  Dimension 

Recall  that  the  Laplace  transform  of  a  function  f(t)  and  its  inverse  transform 
are,  respectively, 

/»oo 

f*(s)  =  /  e~st  f  (t)dt,  (3.23) 

Jo 

and 

-i  pa+ioc 

f(t)  =  x- -  /  estr(s)ds ,  (3.24) 

J  a—ioo 

where  a  >  0  is  arbitrary,  but  must  be  greater  than  the  real  parts  of  all  the  singularities 
of  f*(s).  In  this  thesis,  it  is  assumed  that  f(t)  is  a  real- valued  function,  so  that 
Equation  (3.23)  and  Equation  (3.24)  can  be  replaced  by 


Re{r(s)} 


e  a  f{t)  cos(wt)dt, 


(3.25) 


and 

2£at  roc 

f(t)  = -  /  Re{f*(s)}  cos(wt)dw,  (3.26) 

71  Jo 

where  s  =  a  +  iw.  Solving  Equation  (3.26)  is  generally  not  an  easy  task.  Therefore, 
numerical  techniques  for  approximating  f(t)  are  used. 

One  of  the  most  common  techniques  for  approximating  the  integral  of  Equation 
(3.26)  is  to  exploit  the  relationship  between  the  Laplace  transform  and  the  finite 
Fourier  cosine  transform.  The  paper  by  Dubner  and  Abate  [6]  is  widely  regarded  as 
the  seminal  paper  on  this  technique,  and  will  be  summarized  in  this  section. 
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Consider  a  real  function,  h(t)  such  that  h(t)  =  0  for  t  <  0,  that  has  been 
broken  into  sections,  each  of  length  T.  That  is,  there  are  then  an  infinite  number 
of  intervals  {nT,  (n  +  1)T)  for  n  =  0, 1,  2,  •  •  • .  Each  of  these  is  reflected  through 
its  boundary,  constructing  an  infinite  number  of  even  periodic  functions  gn{t ),  each 
with  period  2 T.  Therefore, 


9n{t) 


hit) 
h{2nT  —  t) 


nT  <t<(n  +  1  )T 
(■ n  —  T)T  <  t  <  nT. 


(3.27) 


Each  gn(t)  is  then  rewritten  so  that  the  functions  are  defined  on  (— T,  T).  For 


n  =  0,2,4,  •••, 

9n(t) 

and  for  n  —  1, 3,  5,  •  •  • , 


h(nT  +  t)  :  0  <  t  <  T 

h(nT  —  t )  :  —  T  <t<  0, 


(3.28) 


9n(t) 


h((n  +  1)T  —  t)  :  0  <  t  <  T 
h((n  +  1)T  +  t )  :  —  T  <t<  0, 


(3.29) 


Taking  the  Fourier  representation  of  each  gn{t ),  and  summing  them  yields 


^9n(t)  = 


n= 0 


Mwo) 


+  ^Aiwk) 


cos 


k= 1 


kir 

Yf' 


(3.30) 


where 

f°°  /  l-TT  \ 

A{wk )  =  I  h{t )  cos  i—tj  dt.  (3.31) 

A(wk)  is  a  Fourier  cosine  transform,  but  if  an  attenuation  factor  is  introduced  by 
letting 

h{t)  =  e-°*f{t),  (3.32) 


it  becomes  the  Laplace  transform  of  the  real  function  f(t)  with  transform  variable 
s  =  a  +  i(k7r/T).  Therefore,  A(w k)  =  Re{/*(s)}.  Multiplying  both  sides  of  Equation 
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(3.30)  by  the  attenuation  factor  eat  yields 

^eatgn(t)  =  ^—  ^Re{/*(a)}  +  ^  Re  (/*  (a+  X  cos  .  (3.33) 

n= 0  l  k= 1  ^  ^  '  ' 

Dubner  and  Abate  [6]  conclude  that  for  any  t  such  that  0  <  t  <  T/2,  the 
inverse  Laplace  transform  can  be  approximated  to  any  desired  accuracy  by 

f(t)  -^eatgn(t)  =  2jr  ^Re{/*(a)}  +  ^  Re  j/*  fa  +  X  cos  ^-t  . 

(3.34) 

Examination  of  Equation  (3.34)  reveals  that  the  approximation  formula  is  simply 
an  application  of  the  trapezoidal  rule  [14]  to  Equation  (3.26).  All  one-dimensional 
inversion  results  in  this  thesis  will  be  obtained  using  algorithms  coded  in  MATLAB®, 
based  on  Equation  (3.34). 

There  are,  however,  other  techniques  to  perform  one-dimensional  inversion. 
Abate,  et  al.  [2]  develop  one  such  technique  in  which  the  desired  function  can 
be  approximated  as  a  weighted  sum  of  Laguerre  functions.  With  some  types  of 
problems,  this  method  performs  better  than  the  Fourier  series  method,  but  in  others 
it  performs  worse  [2],  In  fact,  testing  by  Davies  and  Martin  [5]  on  16  test  problems 
showed  that  the  Laguerre  method  gave  poor  results  on  6  of  those  test  problems.  In 
general,  the  Laguerre  method  is  considered  unsuitable  with  two  types  of  problems. 
It  has  difficulties  with  problems  in  which  the  Laguerre  functions  fail  to  converge 
geometrically,  and  problems  in  which  there  is  geometric  convergence  of  the  Laguerre 
functions,  but  in  which  large  t  values  are  considered. 

3-4-2  Numerical  Inversion  in  Two  Dimensions 

Many  types  of  problems  encountered  in  the  real  world  involve  Laplace  trans¬ 
formations  in  two  dimensions.  In  this  case,  Equations  (3.25)  and  (3.26)  extend 
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naturally  to 

/»oo  /»oo 

f*(s1,s2)=  /  /  e-sitl~s^f(t1,t2)dt1dt2,  (3.35) 

Jo  Jo 

and 

j  /*Ci+200  /»C2+^00 

/  e-iti+-2t2/*(Si>S2)d5id52}  (3.36) 

(z7T?)  Jci—ioo  J  C2—ioo 

where  f(ti,t2)  is  a  real-valued  function  of  t\  and  t2,  f(ti,t2)  =  0  for  t\  or  t2  <  0, 
and  \f(ti,t2)\  <  Me“lfl+Q2t2,  where  aq  and  a2  are  real  numbers  and  M  is  a  positive 
constant.  It  is  also  assumed  that  Ci  >  aq  and  c2  >  a2. 

Moorthy  [15]  provides  a  technique  to  approximate  f(ti,t2)  by  applying  the 
technique  of  Dubner  and  Abate  [6]  in  two  dimensions.  Consequently,  f(ti,t2)  can 
be  approximated  by 

fn(t nh)  =  ^5  |l/*(ci.C2)  +  X]  Re  |/‘  c2  +  ~jr'j  j  cos 

-  Im{/-(cI.c2  +  =)}Sin(^)' 

+  E[r*{/*(c,  +  =c2)}  cc(^l) 

-  Im{r(cI  +  =c2)}Sin(^)‘ 

+  ffWe{rUA2  +  ^)}o0J^  +  ^) 

m=  1  n=l  L  1  \  J  )  \  / 

f  ,■*  f  *n7r  imn\  1  f  rmti  mirt2\ 

+  Ro(/*(c1  +  — ,c2-— ))coS(— ) 

-  Im{/*(c1+=c2  +  =)}Sin(^l  +  !^) 

-  Im{r(c1  +  =c2-=)}Sin(^l  +  ^)]}.  (3.37) 

To  use  the  above  approximation,  it  is  necessary  to  select  values  for  the  param¬ 
eters  ci,  c2,  and  N.  Selecting  an  appropriate  value  for  N  helps  to  control  truncation 
error.  Generally,  N  is  found  by  determining  an  N  for  which  the  difference  between 
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/jv+i  and  f n+n/4  is  negligible.  The  variables  c i  and  C2  can  be  assigned  arbitrarily, 
provided  they  meet  the  previously  mentioned  restriction.  However,  for  many  func¬ 
tions,  it  is  not  easy  to  determine  the  values  of  op  and  a2?  making  the  assignment  of 
ci  and  C2  difficult  as  well. 

As  in  one- dimensional  inversion,  there  are  other  alternatives  to  the  Fourier 
method.  One  such  alternative  was  proposed  by  Abate,  et  al.  [3] ,  in  which  the  Laguerre 
method  is  extended  to  multiple  dimensions.  According  to  the  authors,  the  Laguerre 
method  provides  good  results  for  well-behaved  functions  (i.e.,  the  function  and  its 
derivatives  are  sufficiently  smooth),  but  that  the  Fourier  series  method  proved  to 
be  more  robust.  Like  the  Laguerre  method  in  one  dimension,  the  multidimensional 
extension  performs  poorly  in  problems  with  large  t  values,  and  problems  with  non¬ 
geometric  convergence  of  the  Laguerre  functions. 

In  this  chapter,  a  detailed  description  of  Markov  reward  processes  was  pre¬ 
sented.  Several  results  from  a  transient  analysis  of  such  processes  were  examined, 
including  the  concept  of  a  first  passage  time.  Building  on  the  existing  literature 
on  the  subject,  this  chapter  developed  a  simplified,  analytical  expression  for  the 
Laplace-Stieltjes  transform  of  the  cdf  of  the  first  passage  time  in  one  dimension.  It 
remains  to  demonstrate  this  analytical  result  by  applying  it  to  real-world  operational 
problems,  and  comparing  the  results  with  numerical  results  obtained  by  previously 
established  techniques.  That  will  be  the  focus  of  the  next  chapter. 
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4.  Numerical  Examples 

In  this  chapter,  applications  of  Markov  reward  processes  will  be  considered. 
The  examples  demonstrate  the  usefulness  of  the  techniques  presented  in  this  thesis 
and  involve  problems  from  a  variety  of  real-world  scenarios.  In  some  cases,  the 
results  from  one-dimensional  inversion  will  be  directly  compared  to  results  from 
two-dimensional  inversion. 

4-1  Problems  in  Reliability  Theory 

In  this  section,  the  utility  of  Equation  (3.20)  will  be  demonstrated  through  a 
problem  involving  the  propagation  of  a  crack  in  metallic  materials.  The  problem  will 
be  shown  to  be  a  Markov  reward  process  and  the  technique  developed  in  this  thesis 
will  be  used  to  find  numerical  solutions  for  the  desired  distribution.  These  solutions 
will  be  compared  to  solutions  obtained  in  previous  work  by  Kharoufeh  [11],  obtained 
through  two-dimensional  methods,  as  well  as  with  simulation  results. 

Consider  a  new  metallic  component  which  has  just  been  placed  in  service. 
Initially,  there  are  no  cracks  in  the  component,  but  in  time,  normal  wear  and  fatigue 
cause  the  initiation  of  a  crack.  Continued  operation  causes  the  crack  to  grow.  The 
rate  at  which  this  crack  grows  is  completely  determined  by  the  system’s  random 
environment,  which  consists  of  two  distinct  states.  State  1  causes  the  crack  to  grow 
with  rate  rq  and  state  2  causes  it  to  grow  with  rate  r2-  Furthermore,  it  is  assumed 
the  this  random  environment  can  be  characterized  by  a  continuous-time  Markov 
chain,  {Z(t)  :  t  >  0},  which  alternates  between  the  two  states,  S  =  {1,2},  and  has 
infinitesimal  generator  matrix 


a  —a 
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where  b  is  the  rate  at  which  the  environment  transitions  from  state  1  to  state  2,  and 
a  is  the  rate  at  which  the  environment  transitions  from  state  2  to  state  1. 


This  problem  can  be  considered  as  a  Markov  reward  process.  First,  let  X(t) 
denote  the  length  of  the  crack  at  time  t.  This  is  the  reward  (cost)  in  the  system, 
which  was  previously  denoted  R(t)  in  Chapter  3.  The  set  of  rates  at  which  the  crack 
grows,  {r^rq},  represents  the  set  of  reward  rates,  1Z,  from  Section  3.1.  Finally,  let 
Rd  be  the  diagonal  matrix  of  wear  rates,  which  coincides  with  the  V  matrix  from 
before.  In  this  problem,  T(x)  is  the  random  time  required  for  the  crack  to  first  reach 
a  certain  fixed  length  x  G  M+  (e.g.,  the  time  required  for  the  component  to  fail).  It 
is  the  cumulative  distribution  function  of  the  random  variable  T(x)  that  is  desired. 


For  the  specific  problem  solved  by  Kharoufeh  [11],  the  parameter  values  are  as 
follows.  The  wear  rates  are  rq  =  1.0833  and  r2  =  0.250.  The  generator  matrix  values 
are  a  =  b  =  25/3  and  the  crack  length  threshold  is  set  at  x  —  1.0.  Additionally,  it  is 
assumed,  with  probability  1,  the  system  begins  in  state  1.  Therefore,  the  necessary 
matrices  are 


and 


Q  = 


Rd  — 


-25/3  25/3 

25/3  -25/3 

1.0833  0 

0  0.250 

z0=[l  0], 


5 


Next,  Equation  (3.20)  can  be  used  to  calculate  a  solution  in  the  transform 
space.  Using  the  given  matrices,  that  solution  is 


Gx(s)  — 


[1  0]  exp 


1.0833  0 

0  0.250 


-25/3  25/3 

—  g 

1  0 

25/3  -25/3 

0  1 

1 

1 
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Table  4.1  Cumulative  probability  values  for  a  2-state  reliability  problem. 


X 

t 

2-D  Inversion 

1-D  Inversion 

Simulated 

1.0 

1.197 

0.1259 

0.1270 

0.1288 

1.0 

1.288 

0.2373 

0.2379 

0.2430 

1.0 

1.379 

0.3720 

0.3719 

0.3704 

1.0 

1.470 

0.5128 

0.5122 

0.5140 

1.0 

1.562 

0.6437 

0.6442 

0.6484 

1.0 

1.653 

0.7539 

0.7543 

0.7501 

1.0 

1.744 

0.8396 

0.8397 

0.8360 

1.0 

1.835 

0.9010 

0.9009 

0.9002 

1.0 

1.926 

0.9421 

0.9419 

0.9424 

1.0 

2.018 

0.9677 

0.9679 

0.9698 

1.0 

2.109 

0.9830 

0.9830 

0.9829 

1.0 

2.200 

0.9915 

0.9915 

0.9914 

1.0 

2.291 

0.9958 

0.9959 

0.9957 

1.0 

2.382 

0.9982 

0.9981 

0.9980 

1.0 

2.474 

0.9991 

0.9992 

0.9994 

1.0 

2.565 

0.9995 

0.9997 

0.9996 

1.0 

2.656 

0.9999 

0.9999 

0.9996 

1.0 

2.747 

0.9999 

1.0000 

1.0000 

The  one-dimensional  Laplace  transform  inversion  algorithm  of  Abate  and  Whitt  [1] 
was  used  to  compute  the  analytical  cumulative  distribution  values  at  various  points  in 
time.  These  results  are  shown  in  Table  4.1,  along  with  empirical  results  from  Monte- 
Carlo  simulation.  A  comparison  of  the  results  from  the  one-  and  two-dimensional 
methods  shows  that  the  techniques  provide  similar  results,  which  compare  well  with 
the  benchmark  simulation  results.  This  favorable  comparison  is  easily  seen  in  Figure 
4.1  which  plots  the  distribution  for  each  method. 


Problems  in  Vehicular  Traffic  Flow  Theory 

In  this  section,  two  numerical  examples  will  be  given  in  which  a  first  passage 
time  cumulative  distribution  is  obtained  for  a  vehicle  whose  velocity  is  controlled 
by  a  CTMC.  The  first  example  will  involve  a  two-state  CTMC,  and  the  second 
example  will  involve  a  five-state  CTMC.  Both  examples  were  used  in  Kharoufeh’s 
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t 

Figure  4.1  Cumulative  distribution  of  T(x). 


[8]  derivation  of  two-dimensional  results,  and  will  provide  a  opportunity  to  compare 
one-  and  two-dimensional  results. 

4-2.1  Two-State  CTMC 

In  this  example,  a  vehicle  is  travelling  along  a  straight  line,  and  may  assume 
either  of  two  velocities,  V\  =  50  mph,  or  V2  =  30  mph.  The  particular  velocity 
assumed  at  any  given  time  is  controlled  by  a  two-state  environmental  process,  {Z(t)  : 
t  >  0},  having  state  space  S  =  {1,2}.  One  can  think  of  the  two  states  in  terms  of 
conditions  which  would  normally  affect  a  vehicle’s  velocity.  For  example,  perhaps 
state  1  equates  to  a  clear  day,  while  state  2  represents  rain  or  snow.  Therefore, 
when  Z(t)  is  in  state  i,  the  vehicle  has  velocity  Vj,  i  =  1,2.  The  amount  of  time  the 
system  spends  in  state  1  is  exponentially  distributed  with  rate  /?,  and  the  amount 
of  time  spent  in  state  2  is  exponentially  distributed  with  rate  a.  For  this  problem, 
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oi  =  (5  =  500 hr  1 .  Additionally,  it  is  arbitrarily  assumed  that  the  system  starts  in 
state  1  at  time  0,  so  that  P{Z( 0)  =  1}  =  1. 


Clearly  the  problem  reduces  to  the  basic  Markov  reward  process  described  in 
Chapter  3.  In  this  case,  the  accumulated  reward,  R(t),  is  the  total  distance  travelled 
by  the  vehicle  up  to  time  t.  while  the  set  of  reward  rates,  1Z,  is  the  set  of  velocities 
which  the  vehicle  may  assume.  Finally,  the  environmental  process  in  this  example 
coincides  with  the  stochastic  process,  { Z{t )  :  t  >  0},  in  Section  3.1. 


The  first  step  in  obtaining  the  cumulative  distribution  of  R{t)  is  to  identify 
the  necessary  matrices.  The  infinitesimal  generator  matrix  for  the  CTMC  in  this 
example  is 


-p 

P 

-500 

500 

a 

—a 

500 

-500 

The  V  matrix  is 


Ei 

0 

50 

0 

0 

E 

0 

30 

Because  this  problem  seeks  cdf  values  for  values  of  t  in  minutes  rather  than  hours, 
each  of  these  matrices  is  scaled  by  60.  Finally,  because  the  problem  assumes  (with 
probability  1)  that  the  system  starts  in  state  1  at  time  0,  the  initial  distribution  of 
the  governing  process  is 


zo  =  [1  0], 


Equation  (3.20), 


G^s)  =  z0exp(E  1[Q  -  sl]x)  1, 


can  now  be  applied  to  solve  for  the  cumulative  distribution  at  any  value  of  t  for  a 
particular  value  of  x.  For  example,  if  is  desired  to  find  the  cumulative  distribution 
of  the  random  time  it  takes  to  travel  a  total  distance  of  1  mile  (i.e.,  P{T(1)  <  t}), 
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Table  4.2  Cumulative  probability  values  for  a  2-state  transportation  problem. 


x  (mi) 

t  (min) 

1-D  Inversion 

Simulated 

1.0 

1.25 

0.0036 

0.0033 

1.0 

1.30 

0.0243 

0.0242 

1.0 

1.35 

0.0850 

0.0853 

1.0 

1.40 

0.2019 

0.2006 

1.0 

1.45 

0.3691 

0.3682 

1.0 

1.50 

0.5570 

0.5572 

1.0 

1.55 

0.7278 

0.7283 

1.0 

1.60 

0.8557 

0.8559 

1.0 

1.65 

0.9348 

0.9347 

1.0 

1.70 

0.9754 

0.9762 

1.0 

1.75 

0.9925 

0.9927 

1.0 

1.80 

0.9983 

0.9983 

1.0 

1.85 

0.9998 

0.9996 

1.0 

1.90 

1.0000 

1.0000 

the  above  equation  becomes 
Gx(s)  = 

1 
1 


[1  0]  exp 


50  0 

0  30 


-1  r 


-500  500 

500  -500 


1  0 
0  1 


The  above  equation  can  be  solved,  and  subsequently  inverted  using  an  algo¬ 
rithm  coded  into  the  MATLAB®  software.  The  functions  in  the  code  which  calculate 
the  matrix  exponentiation  in  the  above  equation  utilize  the  Pade  approximation  ap¬ 
proach,  which  is  covered  in  great  depth  in  [20]. 

The  numerical  results  for  various  values  of  t  are  given  in  Table  4.2  along  with 
results  from  Monte-Carlo  simulation.  Comparing  the  results  from  one-dimensional 
inversion  with  the  simulation  results  shows  that  the  one-dimensional  method  gives 
acceptable  results  for  this  problem. 
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4-2.2  Five-State  CTMC 

In  this  example,  a  problem  similar  to  that  of  Section  4.2.1  will  be  examined. 
In  this  case,  the  environmental  process,  {Z{t)  :  t  >  0},  is  a  five-state  CTMC  with 
state  space  S  =  {1,  2,  3, 4,  5}.  When  the  environment  process  is  in  state  i,  the  vehicle 
assumes  a  velocity  of  I /)  =  75 /i  for  i  e  S.  As  in  the  previous  problem,  it  is  assumed 
that,  with  probability  1,  the  system  begins  in  state  1  at  time  0. 

For  this  problem,  the  off-diagonal  entries  of  the  infinitesimal  generator  matrix 
were  assigned  randomly,  based  on  a  uniform  distribution  on  the  interval  (200,400) 
hr-1.  The  resulting  generator  matrix  is 


-919.75 

206.91 

264.85 

238.67 

209.32 

223.01 

-971.71 

301.98 

232.73 

213.98 

343.04 

277.78 

-1283.57 

392.72 

270.03 

353.91 

232.27 

213.69 

-1059.47 

259.59 

370.92 

200.89 

216.80 

225.60 

-1014.21 

The  other  required  matrices  are 

75/1  0  0  0  0 

0  75/2  000 

V  =  o  0  75/3  0  0  , 

000  75/4  0 

0  0  0  0  75/5 

and 

zq  —  [1  0  0  0  0], 

As  before,  the  matrices  are  scaled  by  60,  the  distance  is  set  to  x  —  1.0  mile,  and 
P{T(x )  <  t}  was  computed  for  several  values  of  t.  The  results  are  shown  in  Table 
4.3.  Again,  it  is  evident  that  both  the  one-  and  two-dimensional  techniques  provide 
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Tabic  4.3  Cumulative  probability  values  for  a  5-state  transportation  problem. 


x  (mi) 

t  (min) 

2-D  Inversion 

1-D  Inversion 

Simulated 

1.0 

1.25 

0.0786 

0.0806 

0.0777 

1.0 

1.47 

0.3335 

0.3311 

0.3352 

1.0 

1.70 

0.6859 

0.6922 

0.6865 

1.0 

1.92 

0.9141 

0.9144 

0.9136 

1.0 

2.14 

0.9873 

0.9868 

0.9872 

1.0 

2.37 

0.9991 

0.9990 

0.9991 

1.0 

2.59 

1.0000 

0.9999 

0.9999 

1.0 

2.81 

1.0000 

0.9999 

1.0000 

acceptable  results.  However,  the  amount  of  time  saved  by  using  the  one- dimensional 
technique  grows  dramatically  as  the  number  of  states  in  the  CTMC  increases.  There 
is  no  noticeable  difference  in  computation  time  in  the  one-dimensional  inversion  pro¬ 
cess  between  the  two-state  and  five-state  problems,  whereas  with  the  two-dimensional 
inversion  process,  the  time  required  to  find  a  solution  with  5  states  is  considerably 
greater  than  with  2  states.  In  fact,  it  was  noted  that  the  time  required  to  calculate 
the  cumulative  probability  at  one  point  in  time  in  the  5-state  problem  is  nearly  30 
times  greater  than  the  time  required  in  the  two-state  problem. 


4-3  A  Problem  from  the  Theory  of  Queues 

In  this  section,  a  numerical  problem  will  be  considered  involving  a  fluid  queue. 
For  this  particular  problem,  the  rate  of  fluid  flow  into  an  infinite-capacity  buffer 
(queue)  is  stochastic  and  time- variant. 

Consider  a  sewage  treatment  facility  in  which  the  rate  at  which  sewage  arrives 
is  dependent  on  the  state  of  the  system.  These  states  might  represent  various  levels 
of  disrepair  in  the  system  and/or  blockages  in  the  incoming  pipes.  Assume  that 
there  are  20  possible  states  in  this  system.  In  this  problem,  R(t )  denotes  to  the 
cumulative  amount  of  sewage  flow  into  the  queue  (the  treatment  facility)  up  to  time 
t.  The  state  of  the  environment  at  time  t  is  modelled  by  the  CTMC  { Z(t )  :  t  >  0} 
with  state  space  S  —  {1,  2,  •  •  •  ,  20}.  The  rate  at  which  the  sewage  arrives  is  governed 
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Table  4.4  Cumulative  probability  values  for  a  fluid  queueing  problem. 


X 

t 

1-D  diversion 

Simulated 

1.0 

0.040 

0.0091 

0.0087 

1.0 

0.045 

0.0565 

0.0584 

1.0 

0.050 

0.2090 

0.2083 

1.0 

0.055 

0.4875 

0.4894 

1.0 

0.060 

0.7752 

0.7796 

1.0 

0.065 

0.9402 

0.9417 

1.0 

0.070 

0.9912 

0.9910 

1.0 

0.075 

0.9993 

0.9991 

1.0 

0.080 

1.0000 

0.9999 

by  this  CTMC,  and  the  set  of  possible  rates  is  7 Z  —  {iq,  r2,  •  •  •  ,  r20}  such  that  when 
Z(t)  =  i,  the  sewage  arrives  at  rate  rt .  Modelling  the  system  as  a  Markov  reward 
process  allows  the  time-variant  and  stochastic  nature  of  the  arrival  process  to  be 
incorporated  into  the  analysis.  Suppose  an  analyst  is  interested  in  the  amount  of 
time  required  to  accumulate  1  unit  of  total  flow  into  the  facility  without  regard  to 
its  ultimate  destination.  He  or  she  may  desire  to  find  the  cumulative  distribution  of 
T(l),  the  random  time  required  to  first  accumulate  1  unit  of  fluid. 

Again,  the  first  step  in  finding  a  solution  is  to  identify  the  necessary  matri¬ 
ces.  The  Q  matrix  is  the  infinitesimal  generator  matrix  for  { Z(t )  :  t  >  0}.  For 
this  problem,  the  off-diagonal  entries  of  the  generator  matrix  were  uniformly  dis¬ 
tributed  on  the  interval  (100,  300).  When  the  system  is  in  state  i,  the  rate  at  which 
sewage  arrives  is  rt  =  100/i.  Finally,  it  is  assumed,  with  probability  1,  that  the 
system  begins  in  state  1.  Due  to  the  size  of  these  matrices,  they  are  not  displayed 
in  this  thesis;  however  the  cumulative  probability  results,  along  with  Monte-Carlo 
simulation  results,  are  shown  in  Table  4.4. 

The  major  benefit  of  the  one-dimensional  technique  becomes  more  and  more 
apparent  as  the  dimensionality  of  the  Z  process  increases.  Calculating  the  distribu¬ 
tion  results  in  Table  4.4  via  one-dimensional  inversion  required  approximately  two 
seconds  on  a  Pentium  III  1000  MHz  laptop,  whereas  Monte-Carlo  simulation  of  the 
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Table  4.5  Approximate  Computation  Times  for  Numerical  Examples 


States 

1-D  Inversion 

2-D  Inversion 

2 

2  seconds 

15  seconds 

5 

2  seconds 

7  minutes 

20 

2  seconds 

3  hours 

problem  on  the  same  laptop  requires  approximately  20  minutes  and  two-dimensional 
inversion  would  require  about  3  hours.  Obviously,  as  the  problem  size  increases, 
these  time  discrepancies  are  expected  to  be  even  more  pronounced. 

The  examples  in  this  chapter  clearly  demonstrate  the  usefulness  of  having  a 
closed-form  expression  for  the  Laplace-Stieltjes  transform  of  the  cdf  of  the  first  pas¬ 
sage  time  in  one  dimension.  The  numerical  results  from  this  technique  are  essentially 
equivalent  to  those  obtained  from  the  previous  two-dimensional  technique  (based 
on  comparison  with  benchmark  simulation  results).  However,  the  one- dimensional 
method  provides  tremendous  time  savings  which  are  illustrated  in  table  4.5.  The 
order  of  magnitude  of  this  time  savings  grows  as  the  number  of  states  in  the  CTMC 
grows,  making  it  a  particularly  valuable  technique  for  large  problems. 


4-10 


5.  Conclusions  and  Future  Research 

This  thesis  presented  techniques  for  analyzing  the  cumulative  distribution 
function  of  the  first  passage  time  of  a  Markov  reward  process.  Specifically,  an  exact 
distribution  of  first  passage  times  in  one-dimensional  transform  space  was  developed. 
Additionally,  real-world  applications  of  these  types  of  stochastic,  time-varying  sys¬ 
tems  were  examined,  not  only  to  demonstrate  the  utility  of  this  type  of  modelling, 
but  also  to  provide  numerical  examples  on  which  this  technique  could  be  tested. 

The  first  step  in  the  study  was  to  formally  define  a  Markov  reward  process 
(MRP),  as  well  as  the  concept  of  a  first  passage  time.  The  MRP  was  presented 
as  a  system  consisting  of  a  finite-state  Markov  process  that  accumulates  reward 
over  time.  The  rate  at  which  this  system  accumulates  reward  (or  cost)  is  governed 
directly  by  the  state  of  the  Markov  process.  The  first  passage  time  was  defined 
as  the  random  time  required  to  first  accumulate  a  given  reward  level.  It  was  the 
cumulative  distribution  of  this  random  time  that  was  of  particular  interest  to  this 
study.  A  review  of  the  existing  literature  on  the  subject  revealed  that  solutions  for 
this  distribution  are  typically  found  in  two-dimensional  transform  space.  Numerical 
inversion  is  usually  employed  to  produce  approximations  for  cumulative  distribution 
values.  This  thesis  demonstrated  this  technique  by  reviewing  a  two-dimensional 
transform  space  solution  based  on  the  work  of  Kharoufeh  [8] . 

This  thesis  reduced  the  solution  to  one  dimension  by  two  separate  methods. 
The  first  method  was  a  direct  method  using  an  integrating  factor  to  solve  an  ordinary 
differential  equation.  Previous  work  showed  that  the  distribution  satisfies  the  PDE 
system  shown  in  Equation  (3.6).  By  using  the  Laplace  transform,  the  PDE  system 
was  converted  into  an  ODE  system  with  constant  coefficients.  Through  the  use  of  an 
integrating  factor  and  the  initial  conditions,  a  solution  was  found  for  the  cumulative 
distribution  function  in  one-dimensional  transform  space. 
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The  second  method  involved  using  the  two-dimensional  transform  solution 
found  in  the  existing  literature.  It  was  shown  that  the  two-dimensional  solution 
can  be  analytically  inverted  with  respect  to  the  first  transform  variable,  leaving  a 
solution  in  only  one-dimensional  transform  space.  Numerical  inversion  in  a  single 
dimension  was  then  performed  to  produce  approximations  for  the  desired  cumulative 
probability  values. 

Once  these  results  were  developed,  this  thesis  gave  a  brief  overview  of  some  of 
the  techniques  for  one-  and  two-dimensional  numerical  inversion  of  Laplace  trans¬ 
forms.  These  techniques  provide  the  basis  for  the  MATLAB®  algorithms  which 
were  used  to  produce  numerical  results  in  this  work.  Not  only  did  this  overview 
demonstrate  the  performance  of  numerical  inversion,  but  it  also  helped  to  show  the 
advantage  of  inverting  in  a  single  dimension  rather  than  in  two. 

The  next  step  in  the  thesis  was  to  test  the  one-dimensional  solution  through 
numerical  examples.  An  assortment  of  real-world  scenarios  were  chosen  to  illustrate 
the  utility  of  Markov  reward  processes.  It  is  evident  that  these  processes  can  be 
used  to  model  a  wide  variety  of  systems,  and  that  these  systems  may  have  any 
finite  number  of  states.  This  work  specifically  examined  MRP’s  in  reliability  theory, 
vehicular  traffic  flow  theory,  and  queueing  theory.  Once  each  system  was  modelled 
as  a  MRP,  Equation  (3.20)  was  used  to  find  exact  solutions  in  one-dimensional 
transform  space  for  the  first  passage  time  distributions.  These  were  then  numerically 
inverted  using  Abate  and  Whitt’s  [1]  one- dimensional  inversion  algorithm  to  produce 
numerical  approximations  in  the  time  domain.  In  several  cases,  these  results  were 
compared  to  results  garnered  from  the  two-dimensional  technique  as  well  as  with 
Monte-Carlo  simulation  results.  This  comparison  with  the  simulation  clearly  showed 
that  both  the  one-  and  two-dimensional  methods  produce  favorable  results. 

The  major  contribution  of  this  work  was  the  development  of  Theorem  3.1, 
which  provided  a  solution  for  the  distribution  of  first  passage  times  in  one-dimensional 
transform  space.  To  utilize  this  theorem,  it  is  necessary  to  provide  three  inputs:  the 
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infinitesimal  generator  matrix  for  the  governing  CTMC,  the  diagonal  matrix  of  re¬ 
ward  rates,  and  the  state  of  the  CTMC  at  time  t  —  0.  These  are  the  same  three  inputs 
necessary  to  find  solutions  using  the  existing  two-dimensional  technique.  Therefore, 
no  additional  knowledge  of  the  system  is  required  to  reduce  the  dimensionality  of  the 
problem  from  two  to  one.  In  fact,  less  knowledge  of  the  system  is  required  since  two- 
dimensional  numerical  inversion  algorithms  typically  require  the  analyst  to  supply 
numerous  parameters  based  on  the  estimated  time-domain  function.  The  standard 
one-dimensional  numerical  inversion  algorithms  require  no  such  parameters. 

Perhaps  the  greatest  advantage,  however,  to  finding  solutions  in  one-dimensional 
transform  space  rather  than  two,  is  the  significant  improvement  in  computational 
time.  As  the  size  of  the  governing  CTMC  grows,  the  computation  time  in  the  two- 
dimensional  technique  increases  rapidly.  The  one-dimensional  technique  shows  no 
such  increase.  In  the  examples  presented  in  Chapter  4,  there  was  no  noticeable  dif¬ 
ference  between  the  time  required  to  solve  a  2-state  problem  and  the  time  required  to 
solve  the  20-state  problem.  On  the  other  hand,  with  the  existing  method,  the  com¬ 
puter  required  approximately  1  second  to  solve  for  a  single  cumulative  probability 
value  in  the  2-state  problem,  and  would  need  about  25  minutes  to  solve  for  a  single 
cumulative  probability  value  in  the  20-state  problem.  Thus,  Theorem  3.1  proved  to 
be  useful  in  the  presented  sample  problems,  and  would  be  extremely  convenient  for 
real-world  systems  whose  dimensionality  might  exceed  100  states. 

Techniques  for  finding  the  moments  of  first  passage  times  based  on  the  two- 
dimensional  solution  can  be  found  in  the  existing  literature  [8] .  However,  it  is  possible 
that  a  quicker  and  more  accurate  method  may  exist  as  a  result  of  Theorem  3.1. 
Additionally,  the  results  in  this  work,  as  well  as  any  future  analysis  of  moments, 
could  be  applied  to  larger  problems  (in  excess  of  100  states),  utilizing  real-world 
data.  This  would  clearly  illustrate  how  accurately  MRP’s  can  model  real  systems, 
and  also  show  the  enormous  time  savings  associated  with  one-dimensional  analysis. 
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